diff --git a/tools/Paraview macros/etvk.py b/tools/Paraview macros/etvk.py index ccb79ee1..d7f9458e 100644 --- a/tools/Paraview macros/etvk.py +++ b/tools/Paraview macros/etvk.py @@ -19,6 +19,116 @@ import mmap import os from xml.etree import ElementTree as ET +import json + +def display_pmls_new(pmlthick, dx_dy_dz, nx_ny_nz): + """Display PMLs as box sources using PML thickness values. + Only suitable for gprMax > v.4 + + Args: + pmlthick (tuple): PML thickness values for each slab (cells) + dx_dy_dz (tuple): Spatial resolution (m) + nx_ny_dz (tuple): Domain size (cells) + """ + + pml_names = ['x0', 'y0', 'z0', 'xmax', 'ymax', 'zmax'] + pmls = dict.fromkeys(pml_names, None) + + if pmlthick[0] != 0: + x0 = Box(Center=[pmlthick[0] * dx_dy_dz[0] / 2, + nx_ny_nz[1] * dx_dy_dz[1] / 2, + nx_ny_nz[2] * dx_dy_dz[2] / 2], + XLength=pmlthick[0] * dx_dy_dz[0], + YLength=nx_ny_nz[1] * dx_dy_dz[1], + ZLength=nx_ny_nz[2] * dx_dy_dz[2]) + pmls['x0'] = x0 + + if pmlthick[3] != 0: + xmax = Box(Center=[dx_dy_dz[0] * (nx_ny_nz[0] - pmlthick[3] / 2), + nx_ny_nz[1] * dx_dy_dz[1] / 2, + nx_ny_nz[2] * dx_dy_dz[2] / 2], + XLength=pmlthick[3] * dx_dy_dz[0], + YLength=nx_ny_nz[1] * dx_dy_dz[1], + ZLength=nx_ny_nz[2] * dx_dy_dz[2]) + pmls['xmax'] = xmax + + if pmlthick[1] != 0: + y0 = Box(Center=[nx_ny_nz[0] * dx_dy_dz[0] / 2, + pmlthick[1] * dx_dy_dz[1] / 2, + nx_ny_nz[2] * dx_dy_dz[2] / 2], + XLength=nx_ny_nz[0] * dx_dy_dz[0], + YLength=pmlthick[1] * dx_dy_dz[1], + ZLength=nx_ny_nz[2] * dx_dy_dz[2]) + pmls['y0'] = y0 + + if pmlthick[4] != 0: + ymax = Box(Center=[nx_ny_nz[0] * dx_dy_dz[0] / 2, + dx_dy_dz[1] * (nx_ny_nz[1] - pmlthick[4] / 2), + nx_ny_nz[2] * dx_dy_dz[2] / 2], + XLength=nx_ny_nz[0] * dx_dy_dz[0], + YLength=pmlthick[4] * dx_dy_dz[1], + ZLength=nx_ny_nz[2] * dx_dy_dz[2]) + pmls['ymax'] = ymax + + if pmlthick[2] != 0: + z0 = Box(Center=[nx_ny_nz[0] * dx_dy_dz[0] / 2, + nx_ny_nz[1] * dx_dy_dz[1] / 2, + pmlthick[2] * dx_dy_dz[2] / 2], + XLength=nx_ny_nz[0] * dx_dy_dz[0], + YLength=nx_ny_nz[1] * dx_dy_dz[1], + ZLength=pmlthick[2] * dx_dy_dz[2]) + pmls['z0'] = z0 + + if pmlthick[5] != 0: + zmax = Box(Center=[nx_ny_nz[0] * dx_dy_dz[0] / 2, + nx_ny_nz[1] * dx_dy_dz[1] / 2, + dx_dy_dz[2] * (nx_ny_nz[2] - pmlthick[5] / 2)], + XLength=nx_ny_nz[0] * dx_dy_dz[0], + YLength=nx_ny_nz[1] * dx_dy_dz[1], + ZLength=pmlthick[5] * dx_dy_dz[2]) + pmls['zmax'] = zmax + + # Name PML sources and set opacity + tmp = [] + for pml in pmls: + if pmls[pml]: + RenameSource('PML - ' + pml, pmls[pml]) + Hide(pmls[pml], renderview) + tmp.append(pmls[pml]) + + # Create a group of PMLs to switch on/off easily + if tmp: + pml_gp = AppendDatasets(Input=tmp) + RenameSource('PML - All', pml_gp) + pml_view = Show(pml_gp) + pml_view.Opacity = 0.5 + +def load_src_rx(srcs, dl): + # Display sources and receivers as Paraview box sources + for item in srcs: + p = item['position'] + n = item['name'] + src_rx = Box(Center=[p[0] + dl[0]/2, + p[1] + dl[1]/2, + p[2] + dl[2]/2], + XLength=dl[0], YLength=dl[1], ZLength=dl[2]) + RenameSource(n, src_rx) + Show(src_rx) + + +def load_materials(ml): + # create a threshold filter for each material type + for i, m in enumerate(ml): + threshold = Threshold(Input=model) + threshold.ThresholdRange = [i, i] + threshold.Scalars = ['CELLS', 'Material'] + RenameSource(m, threshold) + + # Show data in view, except for free_space + if i != 1: + thresholddisplay = Show(threshold, renderview) + thresholddisplay.ColorArrayName = ['CELLS', 'Material'] + threshold.UpdatePipeline() from paraview.simple import * @@ -36,50 +146,55 @@ Hide(model) # one file only for now -# Get max numID of materials fp = model.FileName[0] # materials string ms = None -# read the material list comment from the file -b = True +# load json embedded in comment of line 3 with open(fp, 'rb') as f: - while(b): - line = f.readline().decode() - if 'pec,free_space' in line: - # materials string - ms = line - b = False + f.readline() + f.readline() + # comments + c = f.readline().decode() -if ms: - # materials list - ml = ms[5:-5].split(',') +# strip comment tags +c = c[5:-5] +# model information +c = json.loads(c) - # create a threshold filter for each material type - for i, m in enumerate(ml): - threshold = Threshold(Input=model) - threshold.ThresholdRange = [i, i] - threshold.Scalars = ['CELLS', 'Material'] - RenameSource(m, threshold) +# discretisation +dl = c['dx_dy_dz'] +# n voxels +nl = c['nx_ny_nz'] - # Show data in view, except for free_space - if i != 1: - thresholddisplay = Show(threshold, renderview) - thresholddisplay.ColorArrayName = ['CELLS', 'Material'] - threshold.UpdatePipeline() +# load material filters +try: + ml = c['Materials'] + load_materials(ml) +except IndexError: + print('No Materials to load') -# create a threshold filter for the pml -threshold = Threshold(Input=model) -threshold.ThresholdRange = [1, 1] -threshold.Scalars = ['CELLS', 'pml'] -RenameSource('pml', threshold) +# load sources +try: + srcs = c['Sources'] + load_src_rx(srcs, dl) +except IndexError: + print('No sources to load}') -# Show data in view, except for free_space -thresholddisplay = Show(threshold, renderview) -thresholddisplay.ColorArrayName = ['CELLS', 'pml'] -thresholddisplay.Opacity = 0.5 -threshold.UpdatePipeline() +# load rxs +try: + rxs = c['Receivers'] + load_src_rx(rxs, dl) +except IndexError: + print('No receivers to load') + +# load pmls +try: + pt = c['PMLthickness'] + display_pmls_new(pt, dl, nl) +except IndexError: + print('No PMLs to load') RenderAllViews()