krober10nd / SeismicMesh

2D/3D serial and parallel triangular mesh generation tool for finite element methods.
https://seismicmesh.readthedocs.io/
GNU General Public License v3.0
126 stars 32 forks source link

constraining shapes leave some global state #102

Closed nschloe closed 3 years ago

nschloe commented 3 years ago

When trying to mesh the L-shape and something else one after another, SM will error out.

MWE:

import numpy
import SeismicMesh

def l_shape():
    h = 0.1
    bbox = (0.0, 1.0, 0.0, 1.0)

    bbox0 = (0.0, 1.0, 0.0, 0.5)
    rect0 = SeismicMesh.Rectangle(bbox0)

    bbox1 = (0.0, 0.5, 0.0, 1.0)
    rect1 = SeismicMesh.Rectangle(bbox1)

    corners = SeismicMesh.geometry.corners

    def union(x):
        return numpy.minimum.reduce([rect0.eval(x), rect1.eval(x)])

    pfix = numpy.vstack((corners(bbox0), corners(bbox1)))

    points, cells = SeismicMesh.generate_mesh(
        bbox=bbox,
        domain=union,
        edge_length=h,
        pfix=pfix,
        verbose=0,
    )

def ball():
    h = 0.1
    bbox = (-1.0, 1.0, -1.0, 1.0, -1.0, 1.0)

    def ball(x):
        return numpy.sqrt(x[:, 0] ** 2 + x[:, 1] ** 2 + x[:, 2] ** 2) - 1.0

    points, cells = SeismicMesh.generate_mesh(
        bbox=bbox, h0=h, domain=ball, edge_length=h, verbose=False
    )
    points, cells = SeismicMesh.sliver_removal(
        points=points, bbox=bbox, h0=h, domain=ball, edge_length=h, verbose=False
    )

l_shape()
ball()
Constraining 8 fixed points..
Constraining 8 fixed points..
Traceback (most recent call last):
  File "y.py", line 47, in <module>
    ball()
  File "y.py", line 38, in ball
    points, cells = SeismicMesh.generate_mesh(
  File "/home/nschloe/.local/lib/python3.8/site-packages/SeismicMesh/generation/mesh_generator.py", line 399, in generate_mesh
    fh, p, extents = _initialize_points(dim, geps, bbox, fh, fd, h0, opts, pfix, comm)
  File "/home/nschloe/.local/lib/python3.8/site-packages/SeismicMesh/generation/mesh_generator.py", line 756, in _initialize_points
    fh, p, extents = _generate_initial_points(
  File "/home/nschloe/.local/lib/python3.8/site-packages/SeismicMesh/generation/mesh_generator.py", line 741, in _generate_initial_points
    p = np.vstack(
  File "<__array_function__ internals>", line 5, in vstack
  File "/home/nschloe/.local/lib/python3.8/site-packages/numpy/core/shape_base.py", line 283, in vstack
    return _nx.concatenate(arrs, 0)
  File "<__array_function__ internals>", line 5, in concatenate
ValueError: all the input array dimensions for the concatenation axis must match exactly, but along dimension 1, the array at index 0 has size 2 and the array at index 1 has size 3
krober10nd commented 3 years ago

This is an odd one for me. I replicated the behavior. Maybe I'm confused something in the Python language or misunderstanding garbage collection.

When the function exits, the garbage should be collected, no? And then upon the next call the default options should be set and then the kwargs should be parsed. Guess that clearly isn't happening.

krober10nd commented 3 years ago

As a solution, I split the one options dictionary into two dictionaries at the top of generation.py. One for generation options and the other for sliver options.

And then I run the gc.collect() upon exiting the generate_mesh and sliver_removal. This seems to clear things up.

yea that's not a good idea after all.

nschloe commented 3 years ago

Without looking at the code, I'm smelling mutable default arguments.

krober10nd commented 3 years ago

Yes, that's good intuition. This was indeed the problem and as a result I moved the dictionaries that define default options inside their respective functions (in the pull request referenced above). This seemingly fixes the errors.

Actually, this helped me catch some gotcha calls on the README examples (where I didn't pass the bbox to the sliver_removal for instance).

krober10nd commented 3 years ago

Okay, the fix for this was merged in and your MWE now actually works. After I finish the other improvement to the API, I'll release a new version.