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

Bossens-Heckbert smoothing #86

Merged
merged 7 commits into from
Sep 6, 2021
Merged

Bossens-Heckbert smoothing #86

merged 7 commits into from
Sep 6, 2021

Conversation

nschloe
Copy link
Collaborator

@nschloe nschloe commented Sep 6, 2021

Fixes #85.

@@ -263,11 +263,11 @@ def distmesh_smoothing(
# Evaluate element sizes at edge midpoints
edge_midpoints = (mesh.points[edges[:, 1]] + mesh.points[edges[:, 0]]) / 2
p = target_edge_size_function(edge_midpoints.T)
desired_lengths = (
target_lengths = (
f_scale * p * np.sqrt(np.dot(edge_lengths, edge_lengths) / np.dot(p, p))
Copy link

@krober10nd krober10nd Sep 6, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

for the scaling of the target_lengths I like to use the ratio of the medians as it represents the inflation factor better in domains with large disparities in mesh size than the ratio of squares.

# Bossen-Heckbert
def _compute_forces(p, t, fh, min_edge_length, L0mult):
    """Compute the forces on each edge based on the sizing function"""
    N = p.shape[0]
    bars = _get_bars(t)
    barvec = p[bars[:, 0]] - p[bars[:, 1]]  # List of bar vectors
    L = np.sqrt((barvec ** 2).sum(1))  # L = Bar lengths
    L[L == 0] = np.finfo(float).eps
    hbars = fh(p[bars].sum(1) / 2)
    # this line is what I'm referring to in my comment above.
    L0 = hbars * L0mult * (np.nanmedian(L) / np.nanmedian(hbars))
    LN = L / L0
    F = (1 - LN ** 4) * np.exp(-(LN ** 4)) / LN
    Fvec = (
        F[:, None] / LN[:, None].dot(np.ones((1, 2))) * barvec
    )  # Bar forces (x,y components)
    Ftot = _dense(
        bars[:, [0] * 2 + [1] * 2],
        np.repeat([list(range(2)) * 2], len(F), axis=0),
        np.hstack((Fvec, -Fvec)),
        shape=(N, 2),
    )
    return Ftot

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Interesting. I'll add a separate issue for this.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this code in SeismicMesh?

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yea, I was doing it offline. Now I just uploaded some quick implementation in there in krober10nd/SeismicMesh#226

Ideally it should go in maybe as an option. I haven't explored 3D yet.

@nschloe nschloe merged commit 9fbee61 into main Sep 6, 2021
@nschloe nschloe deleted the bossens branch September 6, 2021 17:43
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

potentially better force function
2 participants