Mesh Tensor Generation

[rowanc1]

From @rowanc1 on April 6, 2016 0:29

@lheagy was saying that the mesh tensor needs a bit more work, would be good to do something like:

[ ( 0.05, {"to":100}), (0.05, {"max":25}, 1.3), (25, {"to":1000}) ]

Thoughts on the API?

Copied from original issue: simpeg/simpeg#283

[rowanc1]

This should solve something like the "padding until" problem:

Padding (from MT Tutorial)

author @sgkang

We will use a tensor mesh, which means we can use non-uniform cells for the padding, that is, expanding the width of the cells with depth. We can get away with this because EM fields and fluxes are diffusive and high frequencies will be attenuated as they move through the conductive earth (eg. in seismic, you wouldn't want to do this as the seismic response is dominated by wave propagaion).

We expand by a factor of 1.3 until we are beyond the desired domain extent. Here, we write a small while loop to figure out how many padding cells we should use.

npad = 1  # start with 1 cell
padding_fact = 1.3  # the amount by which we will expand each cell of the padding

def padding_extent(npad):
    """
    given a number of padding cells, this computes how far the padding extends
    """
    padding_widths = cs*padding_fact**(np.arange(npad) + 1)
    return padding_widths.sum()

# keep adding padding until we are beyond the desired extent
padding_z = padding_extent(npad)
while padding_z < domain_extent:
    npad+=1
    padding_z = padding_extent(npad)
    
print("{:1.0f} padding cells extends {:1.2e}m > {:1.2e}m (2 skin depths)".format(
    npad, padding_extent(npad), domain_extent
))

[prisae]

Thanks @lheagy for pointing me to this issue. I might work on this issue a bit later this year, as I will be using it. I think at the moment the signature is a bit cryptic and restricted.

I think there are two different points. (1) is better documenting the signature, and (2) is to improve/extend it. I deal with (1) in this comment, and will write another one for (2).

From going through all the tutorials and examples in the discretize manual I collected the following examples. Are there any more possibilities?

• 1. Providing integers (equidistant grid)

    # Number of cells
    nx = 10
    ny = 20

    M = discretize.TensorMesh([nx, ny])

• 2. Providing vectors

    # Arrays of cell sizes
    hx = np.array([3, 2, 1, 1, 1, 1, 2, 3])
    hy = np.array([3, 1, 1, 3])

    M = discretize.TensorMesh([hx, hy])

• 3. Providing complex input for stretched grid

    # h1 = (
    #    (cellSize, numPad, [, increaseFactor]),
    #    (cellSize, numCore),
    #    (cellSize, numPad, [, increaseFactor])
    # )
    h1 = [(10, 5, -1.3), (5, 20), (10, 3, 1.3)]

    M = discretize.TensorMesh([h1, h1], x0='CN')

As far as I could see, the three cases can be mixed. E.g.,

    hx = [1, 2, 3, 4, 5]
    ny = 10
    M = discretize.TensorMesh([hx, ny])

I haven't gone through the code in detail yet. But how does it distinguish between Case 2 and Case 3? Does it need a value for x0 to trigger Case 3? (Because Case 3 can also be called with only one padded side, eg h1 = ((cellSize, numPad, [increaseFactor])) in which case it is indistinguishable from a 2- or 3-element array of Case 2).

[rowanc1]

• If h is a int, divide the unit space into h cells of size 1/h. Note that the unit space differs on different meshes (e.g. cylinder!).
• If h is a numpy array, this is your spacing vector.
• If h is a list, loop over each element:
  • If h_i is a number, use it as the spacing
  • If h_i is a tuple/list check the length
    • if length 2: append that many sized cells to the overall list. e.g. (2.5, 5) would add 5 cells of size 2.5 to the overall tensor.
    • if length 3: same as above, but change the size of each cell based on the increase factor, if it is negative, flip the list.

The x0 can be used to center or zero the grid easily. The default is to keep all of the axes starting at zero. '000' and [0, 0, '0'] are the same.

[prisae]

Now to (2), extending it. I guess we should keep the existing just as they are, to not break any existing code/example. We could probably add my previous comment to the documentation for clarification.

I would like to see a possibility where you can define (for each, hx, hy, hz):

• min_offset=None
• max_offset=None
• nr_total_cells=None
• stretch_factor=None
• min_cell_size=None
• nr_constant_cells=1

Besides that the parameter x0, but instead of N or C just as a coordinate:

• x0=(0, 0, 0)

The default is x0=(0, 0, 0), but you could set it to any (x, y, z) point. This will permit you to have the origin anywhere in your model.

Edit: as @rowanc1 pointed out, that is already the case.

Assumptions made:

• Center of the model is:
  • (min-offset+max_offset)/2+x0 if min_offset, max_offset are provided;
  • x0 if min_offset, max_offset are not provided.
• Size of the constant cells corresponds to min_cell_size.

This would permit for a couple of different cases, e.g.:

• min_offset, max_offset not provided, everything else is: straightforward calculation of the cells.
• nr_total_cells not provided: cells are created until the first edge that passes min_offset, max_offset.
• stretch_factor not provided: ideal stretch factor is internally calculated to reach with nr_total_cells to min_offset, max_offset.
• min_cell_size not provided: calculated by fitting the other parameters.
• Everything provided: nr_total_cells is then regarded as a maximum. If xmin_offset, xmax_offset with stretch_factor is reached earlier, then less cells are created than nr_total_cells.

[prisae]

Great, @rowanc1, that should be added to the docs too (or is it somewhere already?) If I get around to it I'll do it, but won't be before late September.

[bluetyson]

I am going through the gravity inversion example...however if you want to use a survey that has extents like say 132, -36, 133, -35...eg somewhere in SA and want to define the xo as (132, -36.. how do you then have Z be negative for gravity? I am happy to add to example text when it is clear to me.