API improvement suggestions

[nschloe]

The user interface of SM is more complicated than it needs to be when it comes to bounding boxes, corners etc. Consider the following code for the L-shape:

  def l_shape(h):
      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,
      )
      return points, cells

Suggestion: Define Rectangle objects like

class Rectangle:
     def __init__(self, bbox):
        self.corners = corners
        self.bbox = # ... compute the bbox

     def eval(self, x):
        return # ...

The unions, intersections etc. can compute their bounding boxes and corners from the bounding boxes and corners of the input arguments, e.g.,

class Union:
    def __init__(self, domains: list):
        self.bbox = [
            max(d.bbox[0] for d in domains),
            # ...
        ]

One can check which corners are at the boundary of the union by evaluating the eval() for every domain.

The above code would then look like

  def l_shape(h):
      rect0 = SeismicMesh.Rectangle(0.0, 1.0, 0.0, 0.5)
      rect1 = SeismicMesh.Rectangle(0.0, 0.5, 0.0, 1.0)
      union = SeismicMesh.Union([rect0, rect1])
      points, cells = SeismicMesh.generate_mesh(domain=union, edge_length=h)
      return points, cells

The arguments to Rectangle could be better organized as [0.0, 1.0], [0.0, 0.5].

Defining a new domain would then always require writing a class with eval, a bbox, and (perhaps optionally) corners.

What do you think?

[krober10nd]

Neat that's logical. Yea, corners could be optional as naturally some some geometries don't have them.

I'll get to work sometime today.

[krober10nd]

Okay made some progress over on #106 ....will add the intersection class in a bit. One thing that I'd like to handle is if you want to compute the intersection between say a disk (dim=2) and a cube (dim=3), it could work it out and produce a cylinder.

Maybe the latter isn't a great idea, not sure.

[nschloe]

intersection between say a disk (dim=2) and a cube (dim=3)

Doesn't seem to make sense when first reading it.

I'd say better get union, intersection, and difference down. Perhaps a few primitive shapes like cube, ball, wedge, I don't know. With this, one should be able to do quite a lot already.

[krober10nd]

Yea, I agree. I think that's more of what you would call an extrusion.

[nschloe]

Extrusion, yes, also a good idea. Gmsh has extrude, revolve, twist. For all of these, it's possible to create a 3D model from a 2D one programmatically. That's for the future perhaps.

[jorgensd]

Another minor improvement would be to have an additional keyword argument to plot_sizing_function, such that one can visualize it when using containers etc.

def plot_sizing_function(cell_size, stride=1, comm=None, filename=None):
    # ...
    if filename is None:
        plt.show()
    else:
        plt.savefig(filename)

[krober10nd]

ah yea, that's a good suggestion.

[nschloe]

The following still fails

import SeismicMesh

points, cells = SeismicMesh.generate_mesh(
    h0=0.1,
    domain=SeismicMesh.geometry.Rectangle([-1.0, +1.0, -1.0, +1.0]),
    edge_length=0.1,
)

ValueError: `bbox` must be a tuple

The idea here would be that you don't need to specify a bbox explicitly; it should be taken from domain.

[krober10nd]

The idea here would be that you don't need to specify a bbox explicitly; it should be taken from domain.

It is taken from domain.

https://github.com/krober10nd/SeismicMesh/blob/d71ccee2591b295725afb0abd94c4913fff0b92f/SeismicMesh/generation/mesh_generator.py#L136

It's failing because bbox isn't a tuple

[nschloe]

A tuple as opposed to a list? Why so strict? (Genuine question, perhaps it does make sense.)

[krober10nd]

This should work

 import SeismicMesh

 points, cells = SeismicMesh.generate_mesh(
     domain=SeismicMesh.geometry.Rectangle((-1.0, +1.0, -1.0, +1.0)),
     edge_length=0.1,
 )

[krober10nd]

A tuple as opposed to a list? Why so strict? (Genuine question, perhaps it does make sense.)

Because tuple is immutable and the idea is that the bbox shouldn't mutate.

[nschloe]

Can perhaps also be closed.