Modifying grid.py (appendTriangleBoundary) to different y-max

[Coastal0]

Hello,

I am trying to modify the grid.py file to allow different y-max values on each side of the mesh.

i.e. the left side at sea level (y = 0), while the right side is at topo max (y = ymax)

Currently the yMax is treated as the upper-most node throughout the boundary, however in BERT meshes this primary mesh appears to accurately reflect topography.

Is there a way to replicate this behaviour in pyGimli?

(Originally from https://gitlab.com/resistivity-net/bert/issues/66)

[halbmy]

Hello, I just tried with the simple script

import numpy as np
import matplotlib.pyplot as plt
import pygimli as pg

ex = np.arange(11.)
tx = [ex[0], ex[-1]]
tz = [2., 5.]
ez = np.interp(ex, tx, tz)
plt.scatter(ex, ez)
sensors = [pg.RVector3(exi, ez[i]) for i, exi in enumerate(ex)]
mesh = pg.meshtools.createParaMesh(sensors, boundary=1)
pg.show(mesh)
big = pg.meshtools.appendTriangleBoundary(mesh)

showing the following mesh to be appended:

and observed an error going back to the fact that the bottom edges are detected by the z coordinate being equal to mesh.ymin() which is not the case here. Similar happens to the top edges where only the 7 ones at the right are detected:

bo = [b for b in mesh.boundaries() if b.center().y() == mesh.ymin()]
top = [b for b in mesh.boundaries() if b.center().y() == mesh.ymax()]

If you can detect bottom edges of arbitrary meshes robustly, you might be able to take the next step.

[halbmy]

Let's skip the bottom edge detection (which seems not to be a problem in https://gitlab.com/resistivity-net/bert/issues/66) for a moment and look at the simpler script

import numpy as np
import matplotlib.pyplot as plt
import pygimli as pg

po = [[0, 2],  [0, -3], [10, -3], [10, 4]]
plc = pg.meshtools.createPolygon(po, isClosed=True)
mesh = pg.meshtools.createMesh(plc, quality=34.5, area=3)
big = pg.meshtools.appendTriangleBoundary(mesh)
pg.show(big)

yielding which resembles your original problem. Hence we still need to identify the top nodes robustly to extrapolate their y values to the newly-created top nodes.

Alternatively, we could pass a topographical list (e.g. the sensors or even enhanced by some points outside of it) that is taken for the extrapolation. So we could get rid of the flat (not extrapolated) topo outside...

[halbmy]

I modified grid.py so that the following script works

points = [[0, 2],  [0, -3], [10, -3], [10, 6]]
plc = pg.meshtools.createPolygon(points, isClosed=True)
mesh = pg.meshtools.createMesh(plc, quality=34.5, area=3)
topo = [[-5, 0, 10, 15], [1, 2, 6, 7]]
big = pg.meshtools.appendTriangleBoundary(mesh, topo=topo)
pg.show(big)

It should not do any harm for the existing cases, but I need to check. Topo (consisting of two iterables for x and z) could represent the sensor list, or only the first and last points or any other. @Coastal0 can you check whether this works for your specific case when replacing grid.py from the one out of grid.zip

[Coastal0]

Thanks Thomas, I appreciate the help.

I've ended up with the below (Fig 1). The effects at quite curious at either grid extreme - to the left the y=0 has become slightly elevated (Fig 1-1).

My topography is limited to the parameter mesh (the finer one?), but I will try to add something to extend the topo-list out to the x-min and x-max to see if that solves the issue.

For posterity, the code I've used is below, alongside the concave hull file. The topography finder is a mess, but I think it works.

I had to use matlab to generate the concave hull as I couldn't suss it out quickly enough in python...([http://au.mathworks.com/help/matlab/ref/boundary.html]),

import pygimli as pg
import numpy as np
import os

os.chdir(r"F:\results\QR v13\RemotePC")
fName = 'concaveHullBoundaryVertices.dat' # From MATLAB output.

# %% Topography finder
bVerts = np.loadtxt(fName, dtype=float, delimiter=',',skiprows=1);
bVerts_topo = bVerts[(bVerts[:,1] >= 0)] # get points higher than 0
bVerts_topo = bVerts_topo[np.argsort(bVerts_topo[:,0])]; # sort from 0+

dx = np.lexsort((bVerts_topo[:,1],bVerts_topo[:,0])) # sort by x then z.
u , idx, idc = np.unique(bVerts_topo[dx,0],
                         return_index=1,
                         return_counts = 1,
                         axis=0)
bVtu = bVerts_topo[idx];
for i in enumerate(idc):
    if i[1] > 1:
        print('Found duplicate values')
        bVtu[i[0]] = np.max(bVerts_topo[bVerts_topo[i[0],0] == bVerts_topo[:,0]],0)

# %% Mesh Boundary
FePoly = pg.meshtools.createPolygon(bVerts, 
                                    isClosed = True, 
                                    marker = 98, 
                                    area = 1);
#pg.show(FePoly);
Mesh_FePoly = pg.meshtools.createMesh(poly=FePoly,quality=33)
#pg.show(Mesh_FePoly);
BoundMesh_FePoly = pg.meshtools.appendTriangleBoundary(Mesh_FePoly,
                                                       xbound = 10, 
                                                       ybound = 10, 
                                                       marker = 80, 
                                                       quality = 34,
                                                       isSubSurface = 0,
                                                       topo=bVtu)

pg.show(BoundMesh_FePoly);

concaveHullBoundaryVertices.zip

[carsten-forty2]

just a side note .. appendTriangleBoundary is not supposed for this kind of problems.

appendTriangleBoundary was designed to add a boundary to a grid (a regular mesh of quadrangles) without topography. It might work for an unstructured mesh with topography after somehow hackish work but using a PLC reflecting your whole geometry is the much better approach.

There is already an issue for this to create a solution using the pg.meshtools. In the meantime you can simple use the Bert v.1 tools to create a PLC or the whole mesh.

[Coastal0]

Thanks Carsten.

I haven't been able to work out the Bert tools yet. The pygimli interface is more intuitive for learning the ropes, but I'll keep long looking it.

Ideally I'd like to use the flow as solute modelling similar to Florian's work, and if I can help by contributing I'll gladly do so.

[halbmy]

i agree with Carsten that there is no easy general solution but the changed function works nicely if you transpose your topo and delete the erroneous last point i.e. use topo=bVtu[:-1, :].T

[Coastal0]

Ah wow. So I had the array in the wrong direction...

Thank you greatly for fixing that.

Why was the last point erroneous ...?

[halbmy]

plotting the bVtu points using

plt.scatter(bVtu[:, 0], bVtu[:, 1])

yields which is why I used topo=bVtu[:-1, :].T. I now added a check of the size of topo and transpose it so that it can be both ways. To me it seems not to hurt the original purpose and adds no complexity except one np.interp line which is why I will commit after introducing an exception if the bottom list is not found due to its not straight.

Though I see that in your case it is easier to create the mesh directly using createParaMesh.

[Coastal0]

createParaMesh was indeed easier. I just used the topo vector as the sensor positions, and got a nice result.

I guess I could merge the two meshes together too, somehow!

paramesh = pg.meshtools.createParaMesh(sensors = topo.T)

[carsten-forty2]

which two meshes do you want to merge? your paramesh is all you need for inversion.

[Coastal0]

Once I've worked out how to interpolate my mass-concentration derived resistivity points onto a mesh, I'd like to merge that mesh with the mesh being used for generating the forward model.

My mass-concentration derived point mesh is more-dense than the one generated by the ParaMesh, so I thought I'd need to refine/generate the primary mesh or use the mass-conc. mesh.

[Coastal0]

What was the correct way to merge these meshes.... I'm trying a very basic example and they nodes are just overlying each other!?

(This is trying to setup the basic 2D ERI modelling example in #102 )

[halbmy]

Of course, merging two meshes is not a trivial task, on the contrary it is meaningless to merge meshes that are overlapping each other (how should this work?). You should merge the PLCs (which are also of pygimli mesh type but trivial enough) before meshing.