Cannot see the Complex geometry view using Paraview

[Shikhar Sharma]

Hello Everyone,

I would like to thanks Dr. Craig Warren and Dr. Antonis Giannopoulos for gprMax software it saves lot of time and cost for data collection.

I am a data scientist and have good understanding of  python, I converted below mine simulation code into my_model.in file  for visualization as given below.

#python:
from gprMax.input_cmd_funcs import *

command('title', 'Simulation with a mine object')

# Define simulation parameters
z_dim = 0.002
resolution = 0.002
tsim = 5e-9
B_scan = False

# Set up the domain and grid resolution
domain = domain(x=1.000, y=0.250, z=0.450)
dx = dx_dy_dz(resolution, resolution, z_dim)

# Set the simulation time window
time_window(tsim)

# Define the half-space material
material(permittivity=3.5, conductivity=0.01, permeability=1.0, magconductivity=0, name='bakelite')
material(permittivity=6.0, conductivity=0.01, permeability=1.0, magconductivity=0, name='rubber')
material(permittivity=3.0, conductivity=0.00001, permeability=2.0, magconductivity=0.01, name='dry_sand')
material(permittivity=2.86, conductivity=0.00048, permeability=1.0, magconductivity=9.75, name='TNT')

# Define the waveform
identifier = waveform('ricker', amplitude=1, frequency=1.5e9, identifier='my_ricker')


if B_scan:
    x_ant = 8e-2
else:
    x_ant = domain.x/2 - 1e-3 # in the middle of the x-axis

tx = hertzian_dipole('x',
                     0.100, 0.110, 0.360 , # minus 4 cm in y-direction
                     identifier)

rx(0.100, 0.140, 0.360)     # 2 cm away in x-direction from tx


if B_scan:
    src_steps(dx=0.006)
    rx_steps(dx=0.006)

b0, b1 = box(0, 0, 0,
             1.000, 0.250, 0.200, # same as domain, minus 4 cm in y-direction
             'dry_sand')

#cylinder: 0.500 0.125 0.147 0.500 0.125 0.097 0.002 pec

cylinder(0.500, 0.125, 0.150,
         0.500, 0.125, 0.147,
         0.056, material='rubber')
cylinder(0.500, 0.125, 0.147,
         0.500, 0.125, 0.094,
         0.056, material='bakelite')
cylinder(0.500, 0.125, 0.147,
         0.500, 0.125, 0.097,
         0.053, material='TNT')
cylinder(0.500, 0.125, 0.147,
         0.500, 0.125, 0.097,
         0.002, material='pec')

# Define outputs, geometry view, and snapshots
geometry_view(0, 0, 0, domain.x, domain.y, z_dim, dx.x, dx.y, dx.z, 'my_mine', 'n')

N = 32
for i in range(1, N+1):
    snapshot(0, 0, 0, domain.x, domain.y, z_dim, dx.x, dx.y, dx.z, i*(tsim/N), 'snapshot' + str(i))

#end_python:


And i generated .VTI file using commands 

 python -m gprMax --geometry-only user_models/my_model.in

I am trying  to visualize it using paraview software  so that i can see mine like the given two https://docs.gprmax.com/en/latest/user_libs_landmines.html


Paraview screenshot

Please help me to visualize Landmine . How actually it looks ? Thanks In advance

[Brad Artman]

In my opinion, Paraview is a terrible experience, with a high barrier to entry: developers gone wild with government funding.

Just stick to Python when you can.

I paste below the python functions I use to deal with model.in, model.out,  and vii files. Sorry about the formatting (I use dark mode in PyCharm).

Note in readHDF, one hint for multiple shots in a single HDF, is to make a new pointer to each shot:

src1HDF = srcsHDF['src1'] # make a smaller HDF file from src1 subdictionary

vtiRead gets most of everything from pv.read(file). The pv library can instantly plot and show you all the instances in the file.

import gc

import numpy as np
import pandas as pd
import pyvista as pv
import csv
import h5py
import os
import sys
import warnings
import nexusformat.nexus as nx

def convertLineToArrayOfFloats(s):
for idx, x in enumerate(s):
try:
s[idx] = float(x)
except ValueError:
pass
return np.array(s)
def makeAxis(nSamples, discretization, origin=0):
axis = np.linspace(discretization, (int(nSamples) - 1) * discretization, num=nSamples) + origin
return axis #function makeAxis
def findIndex(location, sampling, origin=0.):
distance = location - origin
n = (distance / sampling) # -1 (minus required?)
return int(n)
def parsemodelin(directory, modelin):
if not os.path.isfile(modelin):
raise Exception("File: " + modelin + " does not exist")
pipe = np.zeros(7)
pipeRadius = 0.

freespace = [1., 0., 1., 0.]
d_freespace = freespace[0]
c_freespace = freespace[1]
nExtenalProperties = np.array([0])
externalProperties = []
with open(modelin, "r") as in_file:
csvreader = csv.reader(in_file, delimiter=' ')
# next(csvreader) # skip top line, but why?
## very poor way to do it, but I need reference to both firs and last lines of file to make below work
for row in csvreader:
if "#title:" in row: #odd I can't make this work
fullPathTitle = row[1]
dirTitle = row[1].split(sep='/')[-1]
dirConcepts = np.array(dirTitle.split(sep='_'))
if row[0] == "#dx_dy_dz:":
samplings = convertLineToArrayOfFloats(row[1:4])
if row[0] == "#rx_array:":
rx_array = convertLineToArrayOfFloats(row[1:])
if "#domain:" in row: #hmm. another way to do it
domains = convertLineToArrayOfFloats(row[1:4])
if row[0] == "#hertzian_dipole:":
source1type = str(row[1])
source1location = convertLineToArrayOfFloats(row[2:5])
if row[0] == "#src_steps:":
sourceSteps = convertLineToArrayOfFloats(row[1:])
if row[0] == "#cylinder:": #see http://docs.gprmax.com/en/latest/input.html?highlight=cylinder#cylinder
pipe = convertLineToArrayOfFloats(row[1:-3])
pipeRadius = pipe[6]
if row[0] == "#box:":
if row[7] == 'lower':
dimensionLower = convertLineToArrayOfFloats(row[1:7])
if row[7] == 'upper':
dimensionUpper = convertLineToArrayOfFloats(row[1:7])
tvd_earth = dimensionUpper[4]
airthickness = (domains[1] - tvd_earth)
airsamples = int(airthickness / samplings[1])
# if homogeneous, default the trench stuff to upper stuff
dimensionBedding = dimensionUpper
dimensionInitial = dimensionUpper
dimensionFill = dimensionUpper
if row[7] == 'bedding':
dimensionBedding = convertLineToArrayOfFloats(row[1:7])
if row[7] == 'initial':
dimensionInitial = convertLineToArrayOfFloats(row[1:7])
if row[7] == 'fill':
dimensionFill = convertLineToArrayOfFloats(row[1:7])
if row[0] == "#material:":
if row[5] == "lower":
materialLower = convertLineToArrayOfFloats(row[1:5])
d_lower = materialLower[0]
c_lower = materialLower[1]
if row[5] == "upper":
materialUpper = convertLineToArrayOfFloats(row[1:5])
d_upper = materialUpper[0]
c_upper = materialUpper[1]
# if homogeneous, default the trench stuff to upper stuff
# not needed as I always included all materials, just didn't use them to model
# materialFill = materialUpper; materialInitial = materialUpper; materialBedding = materialUpper
if row[5] == "bedding":
materialBedding = convertLineToArrayOfFloats(row[1:5])
d_bedding = materialBedding[0]
c_bedding = materialBedding[1]
if row[5] == "initial":
materialInitial = convertLineToArrayOfFloats(row[1:5])
d_initial = materialInitial[0]
c_initial = materialInitial[1]
if row[5] == "fill":
materialFill = convertLineToArrayOfFloats(row[1:5])
d_fill = materialFill[0]
c_fill = materialFill[1]

if row[0] == "seismiclower:":
s_l = (convertLineToArrayOfFloats(row[1:-1]))
externalProperties.append(row[0])
if row[0] == "seismicupper:":
s_u = (convertLineToArrayOfFloats(row[1:-1]))
externalProperties.append(row[0])
if row[0] == "seismicbedding:":
s_b = (convertLineToArrayOfFloats(row[1:-1]))
externalProperties.append(row[0])
if row[0] == "seismicinitial:":
s_i = (convertLineToArrayOfFloats(row[1:-1]))
externalProperties.append(row[0])
if row[0] == "seismicfinal:":
s_f = (convertLineToArrayOfFloats(row[1:-1]))
externalProperties.append(row[0])
if row[0] == "seismicfreespace:":
s_fs = (convertLineToArrayOfFloats(row[1:-1]))
externalProperties.append(row[0])
modelIndictionary = {directory: {"Sampling": samplings, "Domain": domains, "Source Type": source1type,
"Source1 Location": source1location, "Source Steps": sourceSteps, "Source Type": source1type,
"Receiver Line": rx_array, "Domain Extents": domains, "Domain Sampling": samplings,
"Dimension Lower": dimensionLower, "Dimension Upper": dimensionUpper, "Dimension Bedding": dimensionBedding,
"Dimension Fill": dimensionFill, "Dimension Initial": dimensionInitial,
"Material Lower": materialLower, "Material Upper": materialUpper, "Material Bedding": materialBedding,
"Material Initial": materialInitial, "Material Fill": materialFill,
"Dielectric Lower": d_lower, "Conductivity Lower": c_lower,
"Dielectric Upper": d_upper, "Conductivity Upper": c_upper,
"Dielectric Bedding": d_bedding, "Conductivity Bedding": c_bedding,
"Dielectric Initial": d_upper, "Conductivity Initial": c_initial,
"Dielectric Fill": d_fill, "Conductivity Fill": c_fill,
"Dielectric FreeSpace": d_freespace, "Conductivity FreeSpace": c_freespace,
"External Properties": externalProperties, "N External Properies": len(externalProperties),
"Vp Lower": s_l[0], "Vs Lower": s_l[1], "Density Lower": s_l[2],
"Vp Upper": s_u[0], "Vs Upper": s_u[1], "Density Upper": s_u[2],
"Vp Bedding": s_b[0], "Vs Bedding": s_b[1], "Density Bedding": s_b[2],
"Vp Initial": s_i[0], "Vs Initial": s_i[1], "Density Initial": s_i[2],
"Vp Fill": s_f[0], "Vs Fill": s_f[1], "Density Fill": s_f[2],
"Vp FreeSpace": s_fs[0], "Vs FreeSpace": s_fs[1], "Density FreeSpace": s_fs[2],
"Air": np.array((airthickness, airsamples)), "Directory Concepts": dirConcepts,
"Pipe": pipe, "Pipe Radius": pipeRadius,
"TVD_earth": tvd_earth, "Full Path": fullPathTitle, "Directory": dirTitle
}}
gc.collect()
return modelIndictionary
def readVTI(file, getdata=True): #TODO make into dictionary like HDFread above
print("Making model cube from VTI: " + file)
if not os.path.isfile(file):
raise Exception("File: " + file + " does not exist")
vtiRead = pv.read(file)
# vtiRead.plot()
# . . Model dimensions
# WARNING: Jeff's idea of z is depth, whereas gprMax conventions use y for vertical, and z for into plane
# . . Discretization intervals
vSampling = vtiRead.spacing
# . . Define axes
vnx = vtiRead.dimensions[0] - 1 # minus 1 appropriate?
vny = vtiRead.dimensions[1] - 1
vnz = vtiRead.dimensions[2] - 1
vnDimensions = np.array([vnx, vny, vnz])

if getdata:
vtiValues = vtiRead['Material']
modelImage = np.array(vtiValues.reshape((vny, vnx, vnz)))
else:
modelImage = []

#TODO make sure this is the right order of x,y,z
VTIdictionary = {"Sampling": vtiRead.spacing, "NDimensions": vnDimensions,
"Center": vtiRead.center, "Length": vtiRead.length,
"Origin": vtiRead.origin, "Array Names": vtiRead.array_names,
"Another Material Key": vtiRead.active_scalars_name,
"vtiDataCube": modelImage
} #might be some more of interest
gc.collect()
return VTIdictionary # read getVTI
def readHDF(fileHDF, typesHere, getdata=True):
# usefull for exploring HDF tree, but not needed for run time
# filecheck = nx.nxload(fileHDF)
# tree = filecheck.tree
# print(filecheck.tree)
print("Reading-- Headers from: " + fileHDF)
# assumes all coordinate ordinates origins begin at 0.00
fileNameHDF = h5py.File(fileHDF, 'r')
ntHDF = fileNameHDF.attrs['Iterations'] # nt
dtHDF= fileNameHDF.attrs['dt'] # dt
HDFtime = np.hstack((ntHDF, dtHDF))

# domain info
domainSamplingHDF = fileNameHDF.attrs['dx_dy_dz']
domainSizeHDF = fileNameHDF.attrs['nx_ny_nz']
# domain[0]=nx, domain[1]=ny, domain[3]=dx ...

# receiver array info
rnHDF = fileNameHDF.attrs['nrx']
# make a smaller HDF file from rec subdictionary
# is this a pointer or a whole new hdf file? It would be a nice memory trick to not have the data 2x
rxsHDF = fileNameHDF['rxs']
recKeys =rxsHDF.keys()
aa = rxsHDF['rx1']

# read first 2 receiver info to calculate sampling interval
# could probably make references to rather than create new HDF's
rx1HDF = rxsHDF['rx1']
rx2HDF = rxsHDF['rx2']
rxEndHDF = rxsHDF['rx' + str(rnHDF)]
rx1HDFLocationIndexString = rx1HDF.attrs['Name']
rx1ModelLocations = rx1HDF.attrs['Position']
rx2ModelLocations = rx2HDF.attrs['Position']
rxEndModelLocations = rxEndHDF.attrs['Position']
rxdx = (rx2ModelLocations[0] - rx1ModelLocations[0])
rxdy = (rx2ModelLocations[1] - rx1ModelLocations[1])
rxdz = (rx2ModelLocations[2] - rx1ModelLocations[2])
rxDisrectizatin = np.array([rxdx, rxdy, rxdz])
# Assuming regularly sampled linear array in X, no need to query all receiver entries
# Otherwise loop over rx'1' to rx'rnHDF' and make a much bigger out array
# receivers[0]=rec1Origin[m], receivers[1]=rec1Ylocation[m], receivers[3]=rec1Origin[index], receivers[6]=nr
# Don't forget to change floats to integers for Indexes when using
receiversHDF = np.hstack((rx1ModelLocations, rnHDF, rxDisrectizatin, rxEndModelLocations))

# extract source info
nsrcsHDF = fileNameHDF.attrs['nsrc']
srcsHDF = fileNameHDF['srcs'] # make a smaller HDF file from srcs subdictionary
src1HDF = srcsHDF['src1'] # make a smaller HDF file from src1 subdictionary
s1ModelPositionArray = np.hstack((src1HDF.attrs['Position'], nsrcsHDF))
if nsrcsHDF > 1:
src2HDF = srcsHDF['src2']
s2ModelPositionArray = src2HDF.attrs['Position']
sourcesHDF = np.hstack((s1ModelPositionArray, s2ModelPositionArray))
else:
sourcesHDF = s1ModelPositionArray

if sourcesHDF[1] != receiversHDF[1]:
warnings.warn("Source (" + sourcesHDF[1] + ") and Receiver (" + receiversHDF[1] + " not on same Y level")

if getdata:
cube = np.zeros([len(typesHere), rnHDF, ntHDF])
print("Reading-- Data traces from: " + fileHDF)
for itype, typename in enumerate(typesHere):
for iTraceNum in range(rnHDF):
trace = np.array(rxsHDF['/rxs/rx' + str(iTraceNum+1) + '/' + typename][:])
cube[itype, iTraceNum, 0:ntHDF] = trace
else:
cube = []
# if fileNameHDF:
# print("Yes it is open")
# #what are possible ideas?: aaa = fileNameHDF.filename
# if rxsHDF:
# print("rxsHDF exists")
# if srcsHDF:
# print("srcsHDF exists")
fileNameHDF.close()
# if not rxsHDF:
# print("rxsHDF does not exist")
# if not srcsHDF:
# print("srcsHDF does not exist")
# I suppose the child HDF slices go away with the parent
HDFdictionary = {"Domain Sampling": domainSamplingHDF, "Domain Size": domainSizeHDF,
"Receivers": receiversHDF, "Sources": sourcesHDF,
"N receivers": rnHDF, "D receivers": rxdx, "N sources": nsrcsHDF,
"nt": ntHDF, "dt": dtHDF, "Data": cube}
gc.collect()
return HDFdictionary # readHDF
def runStatistics(cube): # asume zero phase, no need to do min and max and delta
print("Finding Big traces in all shots, all types")
nsrcs = cube.shape[0]
ntypes = cube.shape[1]
nrcvs = cube.shape[2]
typeMaxLocation = np.zeros([ntypes, 2], dtype='uint32')
typeMax = np.zeros(ntypes)
for iType in range(ntypes):
for iShot in range(nsrcs):
for iRec in range(nrcvs):
trace = cube[iShot, iType, iRec, :]
#traceMin = min(trace)
traceMax = np.max(trace)
#traceDelta = traceMax - traceMin
# deltascube[iShot, iType, iRec, iRec] = traceDelta
#if traceDelta > typeMax[iType]:
if traceMax > typeMax[iType]:
typeMaxLocation[iType, 0] = iShot
typeMaxLocation[iType, 1] = iRec
#typeMax[iType] = traceDelta
typeMax[iType] = traceMax
# # try max axis
# a = cube[0,0,:,:]
# b = cube[1,0,:,:]
# b = np.max(a, axis=[0, 1], out=buffer)
print(" done Finding")
return typeMaxLocation
def distanceCalc(time, dielectric):
c = 299792458 # speed of light [m/s]
v = c / np.sqrt(dielectric)
distance = v * time # [m]
return distance
def timeCalc(distance, dielectric):
c = 299792458 # speed of light [m/s]
v = c / np.sqrt(dielectric)
time = distance / v # [s]
return time
def semicircle(center_x, center_y, radius, stepsize=0.001):
# ((x - x0) / a) ** 2 + ((y - y0) / b) ** 2 == 1
a = radius
b = radius
x0 = center_x
y0 = center_y
x = np.linspace(-a + x0, a + x0)
y = b * np.sqrt(1 - ((x - x0) / a) ** 2) + y0


return x, y

[Antonis Giannopoulos]

Hi Shikhar,

Thanks for your nice words. The problem is that your geometry view command uses z_dim instead of domain.z So, you are just getting a 2D slice of your model at the lower Z boundary that is the end of the PML and most likely has 0 ID as it should. Also, as you are dealing with small details for targets try using 'f' instead of 'n' to get all the edges in the display and not only volumetric information for the cells. The files will become bigger but you will see more detail like the one I have attached here.

Best 

Antonis

[Brad Artman]

Another frustration in the VTI files is that it's an integer lookup table. So parse the model.in file and apply physical properties with a where() once you figure out the table (like 1 might be free space, 2 one of your box specifications, if you go down the soil model path, it gets quite weird...)

[Craig Warren]

Hi Brad,

Agree with you on Paraview.

Just to answer the lookup table point. This is a design choice with gprMax. There are two approaches you can take when designing an FDTD solver: 

1. One approach is to store the material coefficients at every cell in the FDTD grid. This has the advantage that you have no need for a lookup table, however, it is very memory intensive. Its especially wasteful when you have a relatively limited number of materials in the model, i.e. 100 or so.
2. The other approach is to use a lookup table - which is what we do with gprMax. An integer value is stored for each cell (actually each cell edge - field component) which is used to 'lookup' all the update coefficients for that material. This reduces the memory footprint significantly. Obviously this approach translates to visualising the cell (FDTD mesh geometry). We could convert the mesh to something other than VTK but we have yet to see a viable alternative - which is surprising given how terrible VTK and Paraview are. Always open to listen if there are other geometry visualisation approaches out there for voxelised models.

Kind regards,

Craig