3-d slab implementations using contours

[lhy11009]

@MFraters @danieldouglas92 @alarshi I think this is coming together. I'll use the issue to discuss the key technical steps for the implementation.

The problem is, we are given a slab surface described by depth contours. So each contour is a given depth and composed of a list of points by x and y coordinates. And then we want to query for a mesh in the model domain, for each point in the mesh, figuring out the nearest point on the surface, and the distance of the nearest point from the trench (the contour at 0 depth).

I am basically done with the Python implementation. I don't guarantee this works for every scenario. But it could already handle some sort of complexities.

To move forward, it would be helpful for us to take a look at the steps and their results and provide me with some feedback on:

1. Whether some of the steps could be simplified, or if there is a better approach.
2. Are there any complexities that cannot be handled by these steps?
3. Suggest other test contours or query points?

Also, feel free to take a direct look at the jupyter notebook:
https://github.com/lhy11009/HaMaGeoLib/blob/main/notebooks/experiments/code_implementation/WBContour.ipynb
This could be run by setting the local_dir to one of your system paths

[lhy11009]

Say we are given this contour (a synthetic contour from the sin function):

[lhy11009]

The first step, is to trace from the points trench, with a given angle, to positions on each of the contours. And then trace from the points on each of the contour points backwards to positions on the trench.

a. To trace from the points trench, with a given angle, to positions on each of the contours

Now I assume this angle to be perpendicular to the local trench strike. And by positions, I mean that it's not necessarily the given points, but could be a position, say, on the line segment between two given points.

A straightforward way is as if we are extracting profiles on the map. So here is the results. The points show the mapped points on each of the contours with non-zero depths (5 of them), and the color shows the index of the points at the trench:

The first issue is that there are some overlaps: on each of the contours, the index of points at the trench is not a simple increasing trend. This basically states that occasionally, two points on the trench are mapped to the same location on the contour, which is not surprising given the curve of the trench.

My current fix for this issue, is to enforce that the index of points at the trench must increase and remove the overlaps. Here is the replot of the same figure:

The second issue is that some points on the contours are not mapped to points on the trench. And the contours kind of "shrink" with depth. This issue is actually fixed in the next step.

The third issue is something not shown in this example and one that is a TODO in the implementation. In some other examples, some of the trench points may not have mapped points on the contour. A realistic scenario would be when the trench becomes longer over time.

b. Trace backwards, from points on a given contour, to the positions at the trench. Similarly, by position, it could be on the line segments between given points at the trench.

This is done by performing linear interpolation with the information we collected in the last step, basically:
(index of points at trench, positions on contour) -> (index of points on contour, position of points at trench)

Here is an example with the contour at 600 km depth (deepest one), notice that some of the points don't have mapped trench positions. These points are at the start and end of the contour, related to issue 2 at the last step.

To fix this, I essentially assume the mapped trench could be virtual. As I know where the mapping start to work on the contour and where it doesn't work anymore, I basically take these start and end points, and extract the outside part. These parts are then taken to the surface and become two branches of a virtual trench, which are connected to the start and end points of the trench. Realistically, this would stand for the Cascadia example, where the length of the trench becomes shorter in time.

[lhy11009]

The third step. If we are given a point in the domain, figure out the nearest point on the surface. The white dots in the figure show the test points and their matching nearest points on the surface. This process is split into a few test steps.

a. Build a near-neighbor interpolation with all the points in the contours. With a given test point, figure out the 3 near neighbors in the contours.

b. Decide where the nearest points are on the given contours. Loop the 3 near neighbors, calculate the smallest distance the test point could be connected with another point on the two tine segments on two sides of the near neighbor in their respective contour. Then by comparing the distance we get on every line segment and deciding which line segment actually has the nearest point.

For example, in the figure below, the 3 near neighbors of the first query point (the left most one) are marked by pink colors. And the two dashed pink line segments are where we try to calculate the nearest distance. And the solid pink line segment is actually where we find it in the end.

c. Construct the candidate quads that will have the nearest point on the surface. Another pair of line segments is located on the two adjacent contours to the one we located the pink solid line in the last step. These two are based on where the nearest points are on these two contours.
Then they construct two quads, represented by the red and dark blue quads. The light blue quad serves to illustrate that this method doesn't pre-assume where the quad could be, so if we happened to get an alternative line segment on one of the adjacent contours, we could end up with a different quad.

d. On each of the quad, use triangulation to figure out where the nearest points could be and calculate the distance. Then compare the distance on the two quads and decide where the point could be.

[lhy11009]

And then the third step is to find the related point of the query points on the trench. This step will make use of the information we established in the last two steps, namely the mapping between trench positions and the contour points, as well as the triangle that contains the nearest point to the query points. The same triangulation is performed from the triangle points and their mapped trench positions to the mapped position of the query point.

In the example, the first and last query points are added to show what happens if the trench position is virtual, which stands for a shrinking trench.

[lhy11009]

I have another suggestion regarding the implementation to GWB.

It might be that we can model a lot of code from the current cross-section implementation. And specifically, the contents in this distance_point_from_curved_planes. This would also help if we could store the results in a similar "distance_from_planes" object and pass it to temperature or composition models.

I just took a look at this function when modifying my other PR and realized the layout is a little sophisticated to read. I would like either someone who is familiar with this to explain it to me. And I might try to encapsulate the parts of this function into separate reusable smaller functions for the contour implementation.

[tjhei]

I am not sure I 100% understand what you are trying to do, but given a scalar field, one can compute the gradient, which points into the direction of steepest ascent and is normal to each levelset.
Would that help in determining where "up" is? (note this is then not strictly a line, but a curve)

[lhy11009]

Yes, that is definitely helpful to keep in mind. Especially when we are going to refine the method to path between contours.

I realize that I should add the description of the problem.

The problem is, we are given a slab surface described by depth contours. So each contour is a given depth and composed of a list of points by x and y coordinates. And then we want to query for a mesh in the model domain, for each point in the mesh, figuring out the nearest point on the surface, and the distance of the nearest point from the trench (the contour at 0 depth).