Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Bilinear Remapping #678

Open
aaronzedwick opened this issue Jan 23, 2024 · 4 comments · May be fixed by #1016
Open

Bilinear Remapping #678

aaronzedwick opened this issue Jan 23, 2024 · 4 comments · May be fixed by #1016
Assignees
Labels
new feature New feature or request
Milestone

Comments

@aaronzedwick
Copy link
Member

aaronzedwick commented Jan 23, 2024

Proposed new feature or change:

Implement bilinear remapping in UXarray. Good paper here by David Marsico and Paul Ullrich.

@aaronzedwick aaronzedwick added the new feature New feature or request label Jan 23, 2024
@aaronzedwick aaronzedwick added this to the Remapping milestone Jan 23, 2024
@aaronzedwick aaronzedwick self-assigned this Jan 23, 2024
@aaronzedwick
Copy link
Member Author

function bilinear_remap(source_grid, destination_grid, source_grid_data):
    values = []
    dual = source_grid.get_dual()
    if not source_grid.is_partial():
        # This tree will be used to find the polygon containing the point, as the first polygon that will be searched is one whose face center is nearest
        tree = source_grid.get_ball_tree(coordinates="nodes")

        for point in destination_grid:
            # Find a subset of polygons that contain the point
            polygons_subset = find_polygons_subset(tree, point)

            # Search the subset to find which one contains the point
            for polygon in polygons_subset:
                if point_in_polygon(polygon, point):
                    value = calculate_bilinear_weights(polygon, point)
                    values.append(value)
                    break

    else:  

        tree = source_grid.get_ball_tree(coordinates="nodes")

        for point in destination_grid:
            found = False
            dual_polygons_subset = find_polygons_subset(tree, point)
            
            for polygon in dual_polygons_subset:
                if point_in_polygon(polygon, point):
                    value = calculate_bilinear_weights(polygon, point)
                    values.append(value)
                    found = True
                    break

            # If it isn't found in the dual mesh, search the primal mesh for the point
            if not found:
                primal_tree = source_grid.get_ball_tree(coordinates="face centers")
                primal_polygons_subset = find_polygons_subset(primal_tree, point)
                
                for polygon in primal_polygons_subset:
                    # If found in the primal mesh, assign the nearest neighbor as the value
                    if point_in_polygon(polygon, point):
                        nearest_index = tree.query(point, k=1, return_distance=False)
                        values.append(source_grid_data[nearest_index])
                        found = True
                        break
                        
                # If it is not found in the primal either, assign 0 as the value
                if not found:
                    values.append(0.0)

    return values

function calculate_bilinear_weights(polygon, point):
    if polygon.nodes == 3:
        # Find the weights using barycentric coordinates
        return calculate_barycentric_coordinates(polygon, point)
    elif polygon.nodes == 4:
        # Splitting the quad into two triangles and using barycentric coordinates on the triangle containing the point might be better
        return calculate_weights_with_newton_quadrilaterals(polygon, point)
    else:
        # For polygons that aren't triangles or quads
        sub_triangle = find_sub_triangle(polygon, point)
        return calculate_barycentric_coordinates(sub_triangle, point)

@hongyuchen1030 This is my general idea for the remapping logic.

@aaronzedwick
Copy link
Member Author

aaronzedwick commented Sep 10, 2024

@philipc2 Depending on how we find the subset of polygons we might not use the tree structures at all, right? Except I suppose for the situation where the point is not in a polygon in the dual mesh, but one in the primal. Does our subsetting use tree structures at all?

@hongyuchen1030
Copy link
Contributor

@philipc2 Depending on how we find the subset of polygons we might not use the tree structures at all, right? Except I suppose for the situation where the point is not in a polygon in the dual mesh, but one in the primal. Does our subsetting use tree structures at all?

I am curious for this PR:#938

If this is what I think it is doing, assigning each side of an edge a sign to see which side is within the given face, then we're basically doing a mesh-arrangement. And if we are doing mesh-arrangement, we don't need to do a brute force test about if a point is inside a given face

@aaronzedwick
Copy link
Member Author

@hongyuchen1030 the algorithm will end up using ball tree, so it will be constructed before the call to point in polygon.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
new feature New feature or request
Projects
Status: 🏗 In progress
Development

Successfully merging a pull request may close this issue.

2 participants