3D geometry from DEM and water layer #861
Description
I have a 3D ERT setup with a DEM in the area. And a water layer.
There are 192 electrodes and about 2200 points for topography. All of the included in the input file share.txt.
Water level is 918m. The idea is the following:
I read the file with the coordinates, and I want to create two plc. One that has the topography up to water level (so all points below 918, will be fixed at 918)
data = ert.load('share.txt')
sensors=np.asarray(data.sensors())
sensors[sensors[:,2]<=918,2]=918
And make a plc
plc = mt.createParaMeshPLC3D( sensors,
paraDx=0.33,
paraDepth=30,
paraMaxCellSize=50,
surfaceMeshQuality=34,
marker=1
)
Now I want to make a second plc2, that it will be all inside plc, to avoid intersections. To do that, i remove the xmin,xmax,ymin,ymax
sensors = np.asarray(data.sensors())
ix1=np.where(sensors[:,0]==np.min(sensors[:,0]))[0]
ix2=np.where(sensors[:,0]==np.max(sensors[:,0]))[0]
ix3=np.where(sensors[:,1]==np.min(sensors[:,1]))[0]
ix4=np.where(sensors[:,1]==np.min(sensors[:,1]))[0]
iall=np.unique(np.r_[ix1,ix2,ix3,ix4])
sensors=np.delete(sensors,iall,axis=0)
And now all coordaintes above the water layer, will be shifted by 0.5 meters below, to make sure plc2 is within plc (and make domain extension as small as possible)
sensors[sensors[:,2]>=918,2]=sensors[sensors[:,2]>=918,2]-0.5
plc2 = mt.createParaMeshPLC3D(
sensors,
paraDx=0.33,
paraDepth=30,
paraMaxCellSize=50,
surfaceMeshQuality=34,
paraBoundary=0.1,
boundary=1.5,
marker=5
)
Now merge the two plc
plc=plc+plc2
And make mesh
mesh = mt.createMesh(plc, quality=34)
This gives no error about intersected surfaces, but it reports
Please report this bug to Hang.Si@wias-berlin.de. Include
the message above, your input data set, and the exact
command line you used to run this program, thank you.
FileNotFoundError: [Errno 2] No such file or directory: 'C:\\Users\\marios\\AppData\\Local\\Temp\\tmpl3p219gj-1.node'
Could you adivse how can I solve this (or maybe there is a better way to do it).
IMPORTANT. On the polytools.py, and createParaMeshSurface, the
sZ = pg.interpolate(pntsSurface, pnts[:, 2], surface.positions())
fails returning zeros, making the surface why too complicated so I have added the
sZ[sZ==0]=np.media(pnts[:,2])